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ABSTRACT 

We present a set of three-component Stackel potentials defined by five parameters and 
designed to model the Milky Way. We review the fundamental constraints that any 
model of the Milky Way must satisfy, including the most recent ones derived from 
Hipparcos data, and we study how the parameters of the presented potentials can 
vary in order to match these constraints. Five different valid potentials are presented 
and analyzed in detail: they are designed to be confronted with kincmatical surveys 
in the future, by the construction of three-integral analytic distribution functions. 
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1 INTRODUCTION 

The determination of the mass distribution and dynamical 
structure of the Milky Way is one of the fundamental tasks of 
Galactic Astronomy. Even though the Milky Way is a spiral 
barred galaxy, axisymmetric models are a necessary starting 
point for perturbation analysis and are thus a prerequisite 
if one wants to understand the bar on a theoretical basis. 
K> ; Caldwell & Ostriker (1981), Rohlfs & Kreitschmann (1988) 
» I ' and, recently, Dehnen & Binney (1998) fitted axisymmetric 
mass models of the Milky Way to various measurements of 
the gravitational force field: they concluded that a wide va- 
riety of models can emerge from this fitting process and that 
the mass distribution of the Galaxy is still ill-determined. 

In order to fully exploit kinematical stellar surveys, we 
should construct dynamical models based on Jeans (1915) 
theorem. This theorem states that the phase space distribu- 
tion function of a stellar system in a steady state depends 
only on three isolating integrals of the motion: numerical ex- 
periments (OUongren 1962; Innanen & Papp 1977; Richstone 
1982) showed that most orbits in realistic galactic potentials 
admit three such integrals. The third integral, in addition to 
the binding energy and the vertical component of the angu- 
lar momentum, is not analytic in a general potential: it is 
possible to define an approximate third integral specific to 
particular orbital families (de Zeeuw, Evans & Schwarzschild 
1996; Evans, Hafner & de Zeeuw 1997; De Bruyne, Leeuwin 
& Dejonghe 2000) or to foliate phase space with tori on 
which three action integrals can be defined in a potential 
that difi^ers from the non-integrable one by only a small 
amount (Kaasalainen & Binney 1994; Binney 2002). Instead, 



we choose to construct models with an exact analytic third 
integral by using Stackel potentials (Stackel 1890). These po- 
tentials were introduced into stellar dynamics by Eddington 
(1915) and have since been used in a number of papers (e.g. 
Lynden-Bell 1962; de Zeeuw 1985; Dejonghe 1993; Sevenster, 
Dejonghe & Habing 1995; Durand, Dejonghe & Acker 1996; 
Bienayme 1999) : in fact, the regularity of typical galactic 
potentials can be understood in terms of their proximity to 
Stackel potentials (e.g. Gerhard 1985). 

Our long term goal is to constrain the mass distribu- 
tion of the Galaxy by using kinematical stellar surveys: we 
shall choose a Stackel potential, and use the quadratic pro- 
gramming technique described by Dejonghe (1989) to de- 
termine, for the available surveys, distribution functions in 
the space of the integrals of the motion (see Famaey, Van 
Caelenberg & Dejonghe 2002). Then we shall compare the 
resulting predictions with the surveys. The potential will be 
modified in the light of that comparison. Thus the first step 
is to show that different Stackel potentials are compatible 
with the standard constraints for a mass model of the Milky 
Way. 

In this paper, our goal is to show that a wide variety of 
simple Stackel potentials can fit most known parameters of 
the Milky Way (including Hipparcos latest findings). In or- 
der to do this, we continue the work of Batsleer & Dejonghe 
(1994, hereafter BD) who presented a set of simple Stackel 
potentials with two mass components (halo and disc) and a 
fiat rotation curve, that we generalize by adding a thick disc 
to them since its existence as a separate stellar component 
is now well documented (Ojha et al. 1994; Chen et al. 2001). 
These new potentials are described by five parameters and 
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we will show that many difTerent combinations of these pa- 
rameters are consistent with hmdamental constraints for a 
mass model of the Milky Way. 

In the next section, we review the observational con- 
straints that any mass model of the Milky Way must sat- 
isfy. In section 3, we present the mathematical form of the 
set of Stackel potentials studied in this paper, and we choose 
selection criteria in the light of the constraints reviewed in 
section 2. In section 4, we examine the consequences of those 
selection criteria and we present five potentials differing in 
terms of form and features, and that are plausible for the 
Milky Way. For the conclusions, we refer to section 5. 



2 OBSERVATIONAL CONSTRAINTS 

Recently, Hipparcos data (ESA 1997) have brought nearly 
definitive answers to some long-standing questions in Galac- 
tic Astronomy. In this section, we review the determination 
of those galactic fundamental parameters that any Milky 
Way potential must match. 



2.1 Galactocentric radius of the sun 

The distance of the sun to the galactic center is difficult to 
estimate: the direct method is to compare the average radial 
velocities with the average proper motions of maser spots 
in star forming regions near the galactic center, but this 
method needs very accurate observations and it is affected 
by extinction. If the density of the objects of the stellar halo 
of the Milky Way peaks at the galactic center, the galac- 
tocentric radius of the sun can then also be measured by 
determining the distance of this density peak: the problem 
of this method is that any uncertainty in the absolute magni- 
tude of the stellar candles will reverberate in the estimation 
of the galactocentric distance of the sun. The measurements 
of this distance have been reviewed by Reid (1993), and they 
tend to approach 8 kpc. In this paper, we assume this es- 
timate of 8 kpc for the galactocentic solar radius, with an 
uncertainty of the order of 0.5 kpc. 



2.2 Flat rotation curve 

The stars of the disc travel in nearly circular orbits around 
the galactic center. The determination of the rotation curve 
Wc(cc7), where ru is the galactocentric radius, and in par- 
ticular Wc(w0) is one of the hardest problems in Galactic 
Structure (see section 2.3 for the local shape of the rotation 
curve and the determination of Vciu^Q)). The rotation curve 
is determined by observations of the kinematics of the gas, 
and in particular of the neutral hydrogen 21 cm line, but the 
rotation curve is not well established for the galactic radii 
zu > ujQ . Observations of outer spiral galaxies indicate that 
the rotation curve remains more or less flat after attaining 
its maximum (e.g. Casertano & van Gorkom 1991): the ro- 
tation curve of a mass model of the Milky Way therefore is 
likely to behave similarly. 



2.3 Oort constants and local circular speed 

As the global shape of the rotation curve is not known 
precisely, the determination of its local shape in the solar 



neighbourhood is fundamental for a better knowledge of lo- 
cal Galactic Structure. Lindblad (1925) and Oort (1927a,b) 
developed the model of differential axisymmetric rotation 
with fl = ^ depending only on the distance to the galac- 
tic center. Oort (1927a) introduced two constants (the Oort 
constants A and B), that can be determined from proper 
motions of neighbouring stars and that are directly related 
to the local shape of the rotation curve. Indeed, for a star of 
the solar neighbourhood on a circular orbit, Taylor expand- 
ing Q{vj) to first order in (ro — ^q) yields for the radial 
and transverse velocities (with d denoting the distance to 
the sun): 



Vr = j4dsin2/ 



Vt 

with 

A = 



Adcos2l + Bd 



1 _ dvc 
2^vj Am' 



,B 



'2^zu dm' 



(1) 
(2) 



(3) 



Kuijken & Tremaine (1991) showed that the Taylor ex- 
pansion terms arising from non-circularity of the orbits are 
negligible, just as they should be if the Milky Way is ax- 
isymmetric. Oort (1927b) showed for the first time this si- 
nusoidal effect of the galactic rotation on the radial velocities 
and on the proper motions. He found A = 19kms~^kpc~^ 
and B = — 24kms~^kpc~^. Indeed, by fitting Eq.(^ to ob- 
served proper motions, one can determine A and B if one 
is sure that the frame is not rotating. This last require- 
ment was pretty unsure before the Hipparcos mission. The 
most reliable determination of the Oort constants based 
on the proper motions of the Cepheids that were mea- 
sured by Hipparcos has been derived by Feast & Whitelock 
(1997) who found A = 14.82 ± 0.84km s-^kpc"! and B = 
-12.37 ± 0.64km s^^kpc^V These values were confirmed by 
Mignard (2000) who found A = 14.5 ± 1.0km s"^kpc"^ and 
B = — 11.5±1.0km s^^kpc"^ by using proper motions of dis- 
tant giants. This indicates that the rotation curve is slightly 
declining in the solar neighbourhood. 

The determination of Vc{zuq) follows from the de- 
termination of the Oort constants: using the values of 
Feast & Whitelock (1997), we find Vc{ujq) = (218 ± 
8kms~^)(i:i70/8kpc), a value which is consistent with 
Wc(cc7q) — 220kms~^ that we choose to adopt in this paper. 
Since we impose VcizuQ) — 220kms~^ and zuq — 8±0.5kpc 
for all the potentials defined in section 3, the value of A—B = 



='0 



is in the fixed interval 27.6 ± 1.7km s ^kpc ^ . This 
interval is in accordance with the one deduced from the val- 
ues of A and B determined by Feast & Whitelock (1997), 
i.e. "'^"^9^ = 27.2 ± 0.9kms-^kpc-\ So, the only relevant 
constraint relative to the Oort constants for our potentials 
will be the value of ^^{-cuq). 

However, The Oort constants are not considered as a 
very strong constraint in this paper. Indeed, the measure- 
ment of the proper motion of the compact radio source 
Sgr A* (Backer 1996), seems to indicate that A — B = 
30.1±0.8km s~^kpc~^ when assuming that Sgr A* is station- 
ary with respect to the galactic center. This is inconsistent 
with the determination based on Hipparcos data, and leaves 
us with an uncertainty which still awards a resolution. 
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2.4 Local dynamical mass 

The mass density in the solar neighbourhood is an essen- 
tial constraint for any mass model. It can be surmised that 
this parameter is not directly observed: it is deduced from 
the positions and velocities of tracer stars in the direction 
perpendicular to the galactic plane (the ^-direction) , using 
the Boltzmann and Poisson equations in various forms. 

Jeans equations (first order moments of the Boltzmann 
equation) imply that the vertical acceleration is related 
to the vertical velocity dispersion and to the vertical num- 
ber density g{z) of a population of stars by: 

_2dln{giz)/gi0)) 



dz 



(4) 



if the vertical motion can be separated from the radial and 
azimuthal motion of the stars (Oort-Lindblad approxima- 
tion), if the velocity ellipsoid is aligned with the cylindrical 
coordinate axes (i.e (v^^^Vz) — 0), and if the population is 
isothermal (i.e. is constant as a function of z). 

Oort (1932; 1960) applied this now classical formula to 
late-type stars (A to M) with the assumption that they were 
old enough to have become dynamically well mixed in the 
2-direction. Then he derived using Poisson equation for 
nearly circular orbits: 

47rGp0 = -^-2(A^-B^) (5) 

Oort found p© = 0.09M©pc"^ in 1932 and p© = 
0.15Mqpc~^ in 1960: this last result indicated that there 
might be a lot of dark matter in the disc, by comparison 
with star counts. 

Afterwards, many other similar studies (Yasuda 1961; 
Eelsalu 1961; WooUey & Stewart 1967; Turon Lacarrieu 
1971; Gould & Vandervoort 1972; Jones 1972; Balakirev 
1976; Hill, Hilditch & Barnes 1979) gave discrepant results, 
suffering from inhomogeneities in the data, systematic errors 
due to the use of photometric distances, and undersampling 
near the galactic plane. Because of this undersampling, some 
studies were even based on young O and B stars, assuming 
that the gas and dust out of which they recently formed was 
already roughly relaxed (Stothers & Tech 1964; von Hoerner 
1966). 

More recently, Bahcall (1984a,b,c) described the disc 
matter as a series of isothermal components and analyzed 
the nonlinear self-consistent equations in which the mat- 
ter produces the potential (Poisson equation) and is also 
affected by the potential (Jeans equations for each isother- 
mal component). Bahcall, Flynn & Gould (1992) applied 
this method to a sample of K giants and found p© = 
O.26M0pc~"^, a result leading once again to the presence 
of dark matter in the disc. 

At about the same time, Kuijken & Gilmore (1989a, b,c, 
1991) used another method which does not use the Jeans 
equations, inspired by the method of von Hoerner (1960). 
This method is based on the assumption that the phase 
space distribution function Fz{z,Vz) of any tracer popula- 
tion depends only on 



E.^V-.{z) + -vl 
where 

dv;(z) 



Az 



(6) 



(7) 



This assumption follows from the classical separability of the 
vertical motion of the stars and from the Jeans (1915) the- 
orem in one dimension. The density in configuration space 
g{z) of a tracer population is then related to its density in 
phase space by the integral equation: 



F^{z,v^)Av^ 



/(K) = /(K(^)) 



(8) 



So, there is a unique relation between g{z) and Fz{Ez): Kui- 
jken & Gilmore (1989c, 1991) inverted the Abel transform 
(^) and used the space density profile g{z) of distant K stars 
to predict their velocity distribution at different heights for 
different Vz{z). Then they compared these predictions to the 
velocity data and used a maximum likelihood technique to 
select the best-fitting potential. The data they used were too 
far from the plane to constrain p© and they found a surface 
mass density between z = ±1.1 kpc of 71M0pc~^, a result 
rejecting the presence of dark matter in the disc. 

So, all the investigations between 1932 and the Hippar- 
cos mission had failed to converge to a reliable determination 
of Pq: they were all very uncertain, essentially because of 
inhomogeneities in the tracer samples, undersampling near 
the equatorial plane and the use of photometric distances. 
Hipparcos data solved all these problems: the completeness 
of the Hipparcos survey for stars brighter than m„ — 8.0 
solved the homogeneity problem, the dense probe near the 
plane eliminated the undersampling, while the accurate par- 
allaxes eliminated the use of photometric distances. 

Creze et al. (1998) and Holmberg & Flynn (2000) used 
Eq. (^ to determine pQ from complete samples of nearby 
A-F stars. Given the functions g{z) and Fz{z = 0,v^^q) = 
Fz[Ez) from the observed vertical density and velocity dis- 
tribution, the function Vz(z) can be derived. Then, to deter- 
mine p0 , the Oort constants have to be used in the Poisson 
equation (^: these are also much better known since the 
Hipparcos mission (see section 2.3). Creze et al. (1998) esti- 
mated P0 = 0.076 ±O.O15M0pc"^ while Holmberg & Flynn 
(2000) estimated p© = 0.102 ± O.OlOMgpc"^. These dif- 
ferences between local density estimates using almost the 
same Hipparcos data are due to different a priori hypothesis 
on the vertical potential Vz. Recent investigations at Stras- 
bourg University using a local Hipparcos sample combined 
with two samples at the galactic poles confirm the values 
obtained by Holmberg & Flynn (2000) (Siebert, Bienayme 
& Soubiran 2002). 

We conclude that any mass model of the Milky Way 
must have its solar neighbourhood mass density in the range 
O.OeMepc"^ < P© < O.12M0pc"^ in order to match Hip- 
parcos latest findings (rejecting the presence of a disc-like 
dark matter component). However, it should be remarked 
that all the results presented here are based on the Oort- 
Lindblad approximation, which may not be valid: therefore 
we do not exclude that another subset of potentials, that do 
not match our imposed constraints based on Oort-Lindblad 
approximation, could also be useful for dynamical modeling 
of the Milky Way. 
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3 THE POTENTIALS 

The Stackel potentials are non-rotating potentials for which 
the Hamilton- Jacobi equation is separable and for which all 
orbits admit three analytic integrals of the motion (Stackel 
1890): they form the most general set of potentials that con- 
tain one free function, for which three exact integrals of the 
motion are known, and which can be relevant as models for a 
global potential in galactic dynamics (Lynden-Bell 1962) . In 
this section, we present the mathematical form of the class 
of Stackel potentials that are studied in this paper. 



3.1 Coordinate system 

Axisymmetric Stackel potentials are best expressed in 
spheroidal coordinates {X,(j),iy), with A and u the roots for 
r of 



T + a T + -f 



a < 7 < 0, 



(9) 



and (ccj, (f), z) cylindrical coordinates. The parameters a and 
7 are both constant and we assume them smaller than zero. 
It is convenient to define the axis ratio of the coordinate 
surfaces as e = ^ with a — —a? and 7 — —(? . Together with 
the focal distance A = ^7 — a, the axis ratio defines the 
coordinate system. 



3.2 Three-component Stackel potentials 

An axisymmetric potential is of Stackel form if there exists a 
spheroidal coordinate system (A, 0, in which the potential 
can be written as 



/(A) ^ !(y) 



(10) 



for an arbitrary function /(r) — (t + 7)G(r), G > 0, r = 
A, V. The function — G(A) then represents the potential in 
the z — plane. 

The Milky Way is composed of several mass compo- 
nents: the bulge, the thin disc, the thick disc, the stellar 
halo, the dark halo, and the (dynamically insignificant) in- 
terstellar medium. It is not fundamental that a mass model 
aknowledges explicitly the existence of each of these compo- 
nents. For example, BD presented a set of two-component 
(halo-disc) Stackel potentials with a flat rotation curve. 

Our goal is to show that different Stackel potentials are 
able to fit the latest estimates for the fundamental parame- 
ters of the Galaxy. We first generalize the potentials of BD 
by adding a thick disc to them since its existence as a sep- 
arate stellar component is now well established (Ojha et al. 
1994; Chen et al. 2001). Our potentials have thus three mass 
components: two "flat" components and one spheroidal. The 
spheroidal component accounts for the stellar and dark halo, 
and we shall see in section 4 that our potentials turn out to 
have an effective bulge, which enables us to avoid the explicit 
introduction of a bulge-component. 

We assume that all three components of our potentials 
generate a Stackel potential, with three different coordinate 
systems but the same focal distance. It is straightforward to 
show that the superposition of three Stackel potentials is still 
a Stackel potential when all three coordinate systems have 
the same focal distance. Although the functions /thin, /thick 



and /halo are arbitrary, we assume that they each generate 
a Kuzmin-Kutuzov potential, defined by 

GM 



G{r) = 



(11) 



with M the total mass of the system. Such a potential be- 
comes a point mass potential (V — —^^) in the Galactic 
Plane when A — > 00. We use this potential essentially be- 
cause it is an extremely simple but representative Stackel 
potential: we will show that it is not necessary to use com- 
plicated Stackel potentials in order to match all the known 
and most recently determined parameters of the Milky Way. 

Near the center, in a meridional plane, the lines of con- 
stant mass density corresponding to a Stackel potential are 
approximately ellipsoidal (de Zeeuw, Peletier & Franx 1986). 
For a Kuzmin-Kutuzov potential (see e.g. Dejonghe & de 
Zeeuw 1988), when a > c, the isodensity surfaces are flat- 
tened oblate spheroids, and increasing e = f produces more 
flattening. So, the ratio e has to be high for the thin disc, 
intermediate for the thick disc and close to unity for the 
halo. 

We first define a class of dimensionless potentials Vp in 
dimensionless units {'ujp,Zp), with a focal distance A = 1 
for all three coordinate systems and the central value of 
the potentials equal to —1. Each of these potentials is a 
superposition of three Kuzmin-Kutuzov potentials in three 
different coordinate systems: 

V^(Athin, Athick, Ahalo, t'thin, ^thick, I'halo) = 
/thin(Athin) — /thin (I'thin) 



-fcthin 



Athin ~ I'thin 
/thick (Athick) — /thick (l^thick) 
Athick ^ t'thick 

/halo (Ahalo) — /halo (i-'halo) 



-(1 - fcthin — fcthick) 



A 



halo ^halo 



(12) 



This new class of potentials is thus defined by five param- 
eters (the three axis ratios of the coordinate surfaces and 
the relative contribution of the thin and thick disc masses 
to the total mass, i.e. fcthin and fcthick), which is a reasonable 
augmentation with respect to the BD potentials that were 
defined by three parameters. 
We denote 



Qthin — Qthick ~ 7thin — 7thick = 9l > 
Qthin — Qhalo = 7thin " 7halo = 92 > gi > 0. 



(13) 



So, we can express the class of potentials as a function 
of Athin, fthin and the two constants qi and q2 (we also use 
Eq. (pr])) to give the final form of Vp: 

^^(Athin.J^thin) = -GM( 



vAthin + V i^thin 



+ 



fcthic 



V Athin + (?1 + V^thin + qi 
1 — fcthin — fcthick 



V -^thin + (?2 + V^thin + 92 



(14) 



The dimensionless rotation curve corresponding to such 
potentials is given by: 

vl{-UJp) = TXJpd^pVp{TAjp,Q) 



GM vjp{ 



fcthin 



V Athin (V Athin + Cthin)'^ 
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^thick 



VAthin +9l(\/Athm +91 + \/c?hin +91) 

1 ~ fethin ~ fcthick 



(15) 



\/Athin + I72(\/Athin + + -vA^hiT+ff^)^ 

where -ccip denotes the dimensionless galactocentric radius. 
We shall impose constraints on the shape and flatness of the 
rotation curve in section 3.3. 

In order to transform these dimensionless potentials 
into dimensional ones for the Milky Way, wc denote the di- 
mensionless radius where the rotation curve attains its first 
maximum as Wp^M'- since the Milky Way attains its global 
amplitude of 220 km/s for the first time at a radius of about 
1.5 kpc (Fich & Tremaine 1991), we define a distance scale 
factor rs = kir^rr- The conversion between dimensionless 



and dimensional distances is then given by: 

n7(kpc) = rs'cup 
2 (kpc) = rszp 



(16) 



Then, the total mass M of the Galaxy is adjusted in 
such a way that the dimensional circular velocity at the 
solar radius {ujq = 8 ± 0.5 kpc, see section 2.1) is equal 
to 220 km s^^: wc obtain thus a minimum and a maximum 
value for M, for the two extreme values of the galactocen- 
tric distance of the sun. That adjustment also fixes the local 
mass density in the solar neighbourhood p© (see section 3.3). 

3.3 Selection criteria 

We shall now establish the features that a potential (as de- 
fined in section 3.2) must have to be considered as a plausible 
potential for the Milky Way, in the light of the observational 
constraints reviewed in section 2. 

By definition of the dimensional potentials (sec section 
3.2), the local circular speed Vc{zuq) is equal to 220 kms^^ 
for all the potentials. The first fundamental selection crite- 
rion is the flatness of the rotation curve: this feature can 
be examined in the dimensionless frame-work. Even the BD 
potentials, with only two mass components, could produce 
many different shapes of rotation curves. So, we adopt the 
same simple diagnostic as BD: we denote zUp^M the dimen- 
sionless radius where the circular speed attains its first max- 
imum and we look for a range in where Vc(zu) remains 
larger than 80% of the maximum velocity and is thus more 
or less constant. We denote zup,F the dimensionless radius 
where VcitUp^p) = O-SvciiiUp^M)- A rotation curve is consid- 
ered sufficiently flat if: 

£^ = ^P^Z^i^I^ > 8 (17) 

This is a minimum requirement. 

The second selection criterion is based on the latest de- 
terminations of the Oort constants. For all our potentials. 



^^^^= 27.6 ± 1.7km s-^kpc-^ 



(18) 



which is in accordance with the value determined by Feast 

& Whitclock (1997), ""^"^Q^ = 27.2 ± 0.9km s-^kpc-^ The 
first derivative of the circular velocity in the solar neighbour- 
hood corresponding to the Oort constants found by Feast & 
Whitelock (1997) is fe(ti7©) = -2.4± 1.2km s'^kpc'^ Our 



potentials have two extreme values for this derivative, de- 
pending on the position of the sun; in order to fit the above 
interval, we select the potentials such that: 



max(^(w0)) > -3.6kms"^kpc" 
min(^(ij70)) < -1.2kms-^kpc- 



(19) 



This feature is not as essential as the flatness of the rotation 
curve, because of the intriguing measurement of the proper 
motion of Sgr A* (see section 2.3). 

The last selection criterion is the local mass density in 
the solar neighbourhood: this number is determined by the 
adopted total mass M. We look for its values in both ex- 
treme positions for the sun {zuq, zq) = (7.5, 0.004) kpc and 
{xuq,zq) = (8.5, 0.02) kpc. Following the Hipparcos latest 
flndings reviewed in section 2.4, we select the potentials such 
that: 



max(p0) > O.O6M0PC ^ 
min(p0) < O.12M0pc"^ 



(20) 



This feature is more important than the Oort constants: we 

give the priority to that constraint, contrarily to Sevenster 
et al. (2000) who studied the structure of the inner Galaxy, 
and therefore constructed Stackel potentials with reasonable 
values for the Oort constants but values for p© that are too 
low. 

The total mass of the Galaxy, the mass fractions of the 
discs and the flattening and scale length of the components 
are not well established observationally and are not consid- 
ered as fundamental constraints for the potential. 



4 RESULTS 

Our goal is now to find and select, among the class of po- 
tentials defined in section 3, some representative potentials 
that differ with respect to their form and features, and that 
satisfy the fundamental constraints reviewed in section 2 
and quantified in section 3.3. As it is difficult to visualize 
a five-dimensional parameter space, we shall first visualize 
its structure when ethin and ehaio are fixed (the "winding 
staircase"). Then we shall restrict the parameter space by 
imposing additional constraints on the flatness of the thick 
disc. We shall finally select five representative potentials, 
listed in Table 5. 



4.1 The "winding staircase" 

As an example of the consequences of the choice of the se- 
lection criteria of section 3.3, we look for all the values of 
Ethick, fcthin and fcthick that yield potentials satisfying the se- 
lection criteria for ethin = 75 and ehaio = 1.02 (with the 
thick disc always thicker than the thin disc, i.e. ethick < 75). 
The accordance with the selection criteria results from a 
precise mixing of the two discs and of the halo. Figure 1 
illustrates the volume in parameter space corresponding to 
these satisfactory potentials: the volume looks like a "wind- 
ing staircaise" . When Aithick = 0, we see that all the values 
of ethick are allowed when 13% < fethin < 15%, which results 
from the fact that no thick disc is in fact present. When 
fcthin is close to its maximum possible value (15%), fethick 
has to be zero except when ethick is very close to 1, i.e. when 
the thick disc is a pretty round component (similar to the 
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Figure 1. This "winding staircase" figure displays the possible values of the coordinate axis ratio e = ^ for the thick disc and the 
possible values of the contributions k of the discs to the total mass in order to satisfy the selection criteria of section 3.3 (for fixed values 
ethin = 75 and ehalo = 1-02, and with the thick disc always thicker than the thin disc, i.e. ethick < 75). It gives a rough vision of the 
region of parameter space that satisfies the criteria. Only the region 1.3 < ethick < 2 is in fact relevant for a model of the Milky Way 
(see section 4.2). 



halo) and helps to keep the rotation curve flat. When A;thin 
is smaller than 15%, the possibilities for fcthick are more nu- 
merous, i.e. the thick disc can take a part of the mass. When 
fcthin attains the critical value of 12%, the volume is inflected 
and the thick disc has to be thin and non-zero in order to 
counter-balance the lack of mass in the thin disc. Decreas- 
ing the mass of the thin disc forces the thick disc to become 
thinner and more massive, yielding the "winding staircase" 
volume of Figure 1. The similar volumes in parameter space 
for ethin > 75 have the same form, and are bigger essentially 
because there is more freedom for ethick. 

4.2 Constraints on the scale height of the thick 
disc 

In Figure 1, there are some solutions with a thick disc more 
massive than the thin disc. Even though the mass fractions 
of the discs and the flattening of the components are not 
well established and should be tested in a dynamical study, 
wo know that the mass fraction of the thick disc is smaller 
than that of the thin disc (and represent at most 13% of the 
local thin disc density in the solar neighbourhood) and the 
latest determination of the thick disc scale height based on 
star count data from the Sloan Digital Sky Survey is 665pc 
(Chen et al. 2001). However, Chen et al. (2001) insist on 
the difficulty to converge to a definitive answer: some other 
studies indicate that the scale height could be of the order of 
1 kpc (GUmore 1984; Ojha et al. 1996). We shall reject the 
potentials that are completely inconsistent with those char- 
acteristics (0.6 kpc < hz < 1 kpc) and in particular those 
with a thick disc more massive than the thin disc. 

In order to determine the interval of axis ratios that 
should be considered for the thick disc, we have fitted 



Table 1. Column 1 contains the axis ratio of the coordinate sur- 
faces for the thick disc. Column 2 gives the corresponding scale 
height for a Kuzmin-Kutuzov potential with rs = 1. 
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^z(pc) 


1.3 


914 


1.4 


834 
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782 


1.8 


696 
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663 



the vertical mass distribution corresponding to a simple 
Kuzmin-Kutuzov potential with rs = 1 to an exponential 
law e~^l^^. We conclude that the potentials with 1.3 < 
Ethick < 2 have a scale height between 665 pc and 1 kpc 
and arc the ones we should examine in detail. For each axis 
ratio, Table 1 gives the corresponding scale height. How- 
ever, Figure 2 reveals that the exponential fit is not valid 
for the axis ratios used to model the thin disc, which is not 
surprising since a thin disc could be better understood as 
a superposition of isothermal sheets than by a simple expo- 
nential law. 

4.3 The final selection 

In this section, we look for some three-component potentials 
with different forms and features, all satisfying the selection 
criteria defined in section 3.3 and consistent with what is 
known about the thick disc. 

First of all, we look for the two-component BD po- 
tentials that satisfy the new selection criteria: they are 
listed in Table 2. All the two-component potentials with 
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Table 2. The characteristics of the two-component BD potentials with ejisc = 50,75,130,200 and e^aio = 1.005,1.01,1.02,1.03 have 
been examined. The potentials satisfying the selection criteria are listed in this table (with a step of 0.01 for the relative contribution of 
the disc). Columns 1 and 2 contain the axis ratios of the coordinate surfaces for the disc and the halo. Column 3 contains the relative 
contribution of the disc to the total mass. Column 4 contains the extent of the flat part of the rotation curve (sec Eq. [l7[ ). Column 5 
contains the minimum and maximum local spatial density, while column 6 contains the minimal and maximal local radial derivative of 
the circular velocity, each time for the two extreme positions of the sun. 
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Table 3. The three-component potentials with ethin = 200 and thalo = 1-01 that satisfy the selection criteria are listed in this table 
(with a step of 0.01 for the relative contribution of the two discs). Columns 1, 2, 3 contain the axis ratios of the coordinate surfaces for 
the two discs and the halo. Columns 4 and 5 contain the relative contribution of the repectively thin and thick disc to the total mass. 
Column 6 contains the extent of the flat part of the rotation curve (see Eq. . Column 7 contains the minimum and maximum local 
spatial density, while column 8 contains the minimal and maximal local radial derivative of the circular velocity, each time for the two 
extreme positions of the sun. 
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Figure 2. Profile of the logarithm of vertical density at -n? = 8 kpc for Kuzmin-Kutuzov potentials with = 1 and e = 1.3 (top left), 
e = 2 (top right), e = 75 (bottom left) and e = 200 (bottom right). Only the two first cases resemble exponentials. 



Table 4. The three-component potentials with cthin = 75 and Chalo = 101 that satisfy the selection criteria are listed in this table (with 
a step of 0.01 for the relative contribution of the two discs). Columns have the same meaning as in table 3. 
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Table 5. Among the class of potentials defined in section 3, five different potentials regarding form and features have been selected. 
Columns 1, 2, 3 contain the axis ratios of the coordinate surfaces for the two discs and the halo. Columns 4 and 5 contain the relative 
contribution of the repectively thin and thick disc to the total mass. Column 6 contains the extent of the fiat part of the rotation curve 
(see Eq. [l7| ). Column 7 contains the scale factor which corresponds to the focal distance of the coordinate system of the dimensional 
potential. Column 8 contains the minimum and maximum local spatial density in Mqpc~^, while column 9 contains the minimum and 
maximum local radial derivative of the circular velocity in kms~^kpc~^ and column 10 the minimum and maximum total mass of the 
Galaxy in 10^^ Mq, each time for the two extreme positions of the sun. The scale length in the equatorial plane down to 3 kpc has been 
calculated and is presented in column 11 (in kpc). 
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Figure 3. Mass isodensity curves in a meridional plane for the five potentials of Table 5, with two different scales (left panel: zoom on 
the disc, right panel: large scale view of the halo). 
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Figure 4. The logarithm of the mass density in the equatorial plane for the two potentials with extreme scale lengths. These curves very 
much resemble each other, and the effective bulge appears clearly. 



Cdisc = 50,75, 130,200 and ehaio = 1.005, 1.01, 1.02, 1.03 are 
examined. We do not consider discs with o/c > 200 because 
then the uncertainty on p© becomes too large. We see in 
Table 2 that, in order to reproduce the Oort constants in 
the two-component framework, the shape of the halo can- 
not vary (ehaio = 1-02). 

If wo take the two-component potentials as a starting 
point, there are two ways to add a thick disc. The first way 
is to decrease the contribution of the halo and to put the 
remaining mass into the thick disc: the local density in the 
solar neighbourhood is then slightly larger while the rota- 
tion curve is decreasing faster. If wc take the third potential 
of Table 2 as a starting point, the first potential of Table 
5 (potential I) illustrates this first case. The other way is 
to decrease the contribution of the thin disc and to put the 
remaining mass into the thick disc: the local density is then 
slightly decreasing while the rotation curve is more flat. If 
we take the fourth potential of Table 2 as a starting point, 
the first potential of Table 5 illustrates this second case: this 
potential (potential I) is a three-component potential satis- 
fying the selection criteria and is selected to be analyzed in 
detail and confronted with kinematical surveys in the future. 

The presence of a third component allows more freedom 
for the shape of the halo, so we look for three-component 
potentials with a halo rounder than ehaio = 1-02. In order to 
keep the rotation curve fiat and retain the local density as 



well as the Oort constants in the allowed interval, we need 
to couple a very thin disc with the rounder halo: indeed, our 
investigations show that no solution can be found for ethin = 
50 and ehaio = 1.01. However, if we take ethin = 200, Table 
3 gives solutions for a halo with ehaio = 1.01: we select the 
solution where the mass of the thick disc relative to the thin 
disc is the smallest (potential III of Table 5). Remark that a 
similar Table for ehaio = 1-02 would contain 292 entries and 
is omitted here. There are much less solutions when ethin ~ 
75, as can be seen in Table 4. However, as stated in section 
3.3, we do not assign high priority to the Oort constants, 
and we select a potential with ethin = 75, ehaio = 1-01 and a 
relative mass of the thick disc relative to the thin disc of 13% 
(i.e. a smaller fraction than aiij' tjf tlic solutions of Table 4), 
but with a quite large local radial derivative of the circular 
speed (potential IV of Table 5). 

For Thaio = 1.005, it is totally impossible to find a po- 
tential satisfying the Oort constants criterion: the radial 
derivative of the circular speed in the solar neighbourhood 
is always positive. Nevertheless, if one is willing to ignore 
the estimates of A and B, one could select a potential with 
ehaio = 1.005 and reasonably low values for the cirular speed 
radial derivative in the solax neighbourhood (potential V of 
Table 5). 

Finally, we select a potential (potential II) satisfying 
all the criteria, for which the interval in is precisely 
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Figure 5. The rotation curves of the five selected potentials of Table 5. The total mass used to plot thes curves is the mean total mass 
of the two extreme values of Table 5. 



[0.06Afgpc"^, 0.12MqPc"^], and which is close to the Chen 
et al. (2001) findings , i.e. ethin = 200 which is the thinnest 
thin disc that we consider in order not to have a too large 
interval for pQ, a relative mass of the thick disc to the thin 
disc of 10%, a scale height of the thick disc of 612.5 pc and 
a relatively large Ep (Ep = 13.44). 

Table 5 summarizes the main features of the selected 

Stackel potentials with different forms and features and that 
we shall use for dynamical modeling of the Milky Way: po- 
tentials III and V have a very thin disc associated with 
a quite massive thick disc and a pretty round halo, while 
potentials I and IV have a thicker thin disc with a quasi- 
negligible thick disc (Figure 3 shows the mass isodensity 
curves of each potential in a meridional plane for two differ- 
ent scales). It should be noted that the total masses asso- 
ciated with those potentials are very different and become 
larger with a rounder halo and that a rounder halo implies 
that this halo is much more extended. A closer look to the 
mass density in the equatorial plane indicates that, for each 
potential, the density grows faster than an exponential in 
the central 3 kpc corresponding to the bulge region: the po- 
tentials have thus an effective bulge, which did enable us to 
avoid the introduction of an explicit bulge component. Wo 
have fitted the mass density in the plane to an exponen- 
tial law down to zu = 3 kpc in order to check that the scale 
length of the disc is realistic: the last column of Table 5 gives 



the scale length corresponding to each selected potential and 
we conclude that they are realistic but do not distinguish the 
different potentials. For the potentials with the biggest and 
smallest scale length, Figure 4 illustrates the shape of the 
logarithm of the density in the plane (the other potentials 
have a similar shape for that curve) . Finally, Figure 5 shows 
the rotation curve associated with each of the five selected 
potentials: the rotation curve of potential II is more flat than 
the one of potential I in the vicinity of the sun, while the ro- 
tation curves of the potentials with ehaio = 1.01 (potentials 
III and IV) are even more flat and the one of potential V is 
slightly increasing. 



5 CONCLUSIONS 

In this paper, we have shown that some different simple 
Stackel potentials can fit most known parameters of the 
Milky Way (especially Hipparcos latest findings). First, we 
have reviewed the galactic fundamental parameters that 
any Milky Way potential must match and that investiga- 
tions following the Hipparcos mission have readjusted. Then, 
we have generalized the two-component BD potentials by 
adding a thick disc, and we have studied how the parame- 
ters can vary in order to satisfy selection criteria based on 
the latest observational constraints. We have shown that the 
presence of a thick disc allows more flexibility in the form 
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of the potentials, especially for the shape of the halo and 
we have selected five different valid potentials listed in Ta- 
ble 5. It should be noted that, in fact, there could be two 
different thick discs in the Galaxy, a very thick one (Chiba 
& Beers 2000; Gilmore, Wyse & Norris 2002) and a smaUer 
one (Soubiran, Bienayme & Siebert 2002): in that case, the 
three-component modelling presented in this paper could be 
eeisily extended, but this would imply a growth of parameter 
space. 

A major result of this paper is that, even though Stackel 
potentials are negligible in the set of all potentials, many of 
them are still able to match the most recent estimates for the 
parameters of the Milky Way, and furthermore very simple 
ones (superpositions of three Kuzmin-Kutuzov potentials) 
are sufficient to do this. 

The potentials defined in this paper have already been 
used in Famaey et al. (2002). We plan to further validate 
each of the five proposed potentials by confronting them 
with kinematical surveys. This we shall do by constructing 
three-integral analytic distribution functions in those poten- 
tials, using a quadratic programming technique (Dejonghe 
1989; for an overview see Dejonghe et al. 2001). 
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ABSTRACT 

We present a set of three-component Stackel potentials defined by five parameters and 
designed to model the Milky Way. We review the fundamental constraints that any 
model of the Milky Way must satisfy, including the most recent ones derived from 
Hipparcos data, and we study how the parameters of the presented potentials can 
vary in order to match these constraints. Five different valid potentials are presented 
and analyzed in detail: they are designed to be confronted with kinematical surveys 
in the future, by the construction of three-integral analytic distribution functions. 

Key words: 

Galaxy: kinematics and dynamics - Galaxy: stucture 



1 INTRODUCTION 

The determination of the mass distribution and dynamical 
structure of the Milky Way is one of the fundamental tasks of 
Galactic Astronomy. Even though the Milky Way is a spiral 
barred galaxy, axisymmctric models are a necessary starting 
point for perturbation analysis and are thus a prerequisite 
if one wants to understand the bar on a theoretical basis. 
Caldwell & Ostriker (1981), Rohlfs & Kreitschmann (1988) 
and, recently, Dchncn & Binncy (1998) fitted axisymmctric 
mass models of the Milky Way to various measurements of 
the gravitational force field: they concluded that a wide va- 
riety of models can emerge from this fitting process and that 
the mass distribution of the Galaxy is still ill-determined. 

In order to fully exploit kinematical stellar surveys, we 
should construct dynamical models based on Jeans (1915) 
theorem. This theorem states that the phase space distribu- 
tion function of a stellar system in a steady state depends 
only on three isolating integrals of the motion: numerical ex- 
periments (Ollongren 1962; Innanen & Papp 1977; Richstone 
1982) showed that most orbits in realistic galactic potentials 
admit three such integrals. The third integral, in addition to 
the binding energy and the vertical component of the angu- 
lar momentum, is not analytic in a general potential: it is 
possible to define an approximate third integral specific to 
particular orbital families (de Zeeuw, Evans & Schwarzschild 
1996; Evans, Hafner & de Zeeuw 1997; De Bruyne, Leeuwin 
& Dejonghe 2000) or to foliate phase space with tori on 
which three action integrals can be defined in a potential 
that difi^ers from the non-integrable one by only a small 
amount (Kaasalainen & Binney 1994; Binney 2002). Instead, 



we choose to construct models with an exact analytic third 
integral by using Stackel potentials (Stackel 1890). These po- 
tentials were introduced into stellar dynamics by Eddington 
(1915) and have since been used in a number of papers (e.g. 
Lynden-Bell 1962; de Zeeuw 1985; Dejonghe 1993; Scvcnster, 
Dejonghe & Habing 1995; Durand, Dejonghe & Acker 1996; 
Bienayme 1999) : in fact, the regularity of typical galactic 
potentials can be understood in terms of their proximity to 
Stackel potentials (e.g. Gerhard 1985). 

Our long term goal is to constrain the mass distribu- 
tion of the Galaxy by using kinematical stellar surveys: we 
shall choose a Stackel potential, and use the quadratic pro- 
gramming technique described by Dejonghe (1989) to de- 
termine, for the available surveys, distribution functions in 
the space of the integrals of the motion (see Famaey, Van 
Caelenberg & Dejonghe 2002). Then we shall compare the 
resulting predictions with the surveys. The potential will be 
modified in the light of that comparison. Thus the first step 
is to show that different Stackel potentials are compatible 
with the standard constraints for a mass model of the Milky 
Way. 

In this paper, our goal is to show that a wide variety of 
simple Stackel potentials can fit most known parameters of 
the Milky Way (including Hipparcos latest findings). In or- 
der to do this, we continue the work of Batsleer & Dejonghe 
(1994, hereafter BD) who presented a set of simple Stackel 
potentials with two mass components (halo and disc) and a 
fiat rotation curve, that we generalize by adding a thick disc 
to them since its existence as a separate stellar component 
is now well documented (Ojha et al. 1994; Chen et al. 2001). 
These new potentials are described by five parameters and 
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we will show that many different combinations of these pa- 
rameters are consistent with fundamental constraints for a 
mass model of the Milky Way. 

In the next section, we review the observational con- 
straints that any mass model of the Milky Way must sat- 
isfy. In section 3, we present the mathematical form of the 
set of Stackel potentials studied in this paper, and we choose 
selection criteria in the light of the constraints reviewed in 
section 2. In section 4, we examine the consequences of those 
selection criteria and we present five potentials differing in 
terms of form and features, and that are plausible for the 
Milky Way. For the conclusions, we refer to section 5. 



neighbourhood is fundamental for a better knowledge of lo- 
cal Galactic Structure. Lindblad (1925) and Oort (1927a,b) 
developed the model of differential axisymmetric rotation 
with n = depending only on the distance zu to the galac- 
tic center. Oort (1927a) introduced two constants (the Oort 
constants A and S), that can be determined from proper 
motions of neighbouring stars and that are directly related 
to the local shape of the rotation curve. Indeed, for a star of 
the solar neighbourhood on a circular orbit, Taylor expand- 
ing n('C!7) to first order in {zu — taJq) yields for the radial 
and transverse velocities (with d denoting the distance to 
the sun): 



2 OBSERVATIONAL CONSTRAINTS 

Recently, Hipparcos data (ESA 1997) have brought nearly 
definitive answers to some long-standing questions in Galac- 
tic Astronomy. In this section, we review the determination 
of those galactic fundamental parameters that any Milky 
Way potential must match. 



2.1 Galactocentric radius of the sun 

The distance of the sun to the galactic center is difficult to 
estimate: the direct method is to compare the average radial 
velocities witli the average proper motions of maser spots 
in star forming regions near the galactic center, but this 
method needs very accurate observations and it is affected 
by extinction. If the density of the objects of the stellar halo 
of the Milky Way peaks at the galactic center, the galac- 
tocentric radius of the sun can then also be measured by 
determining the distance of this density peak: the problem 
of this method is that any uncertainty in the absolute magni- 
tude of the stellar candles will reverberate in the estimation 
of the galactocentric distance of the sun. The measurements 
of this distance have been reviewed by Reid (1993), and they 
tend to approach 8 kpc. In this paper, we assume this es- 
timate of 8 kpc for the galactocentic solar radius, with an 
uncertainty of the order of 0.5 kpc. 



2.2 Flat rotation curve 

The stars of the disc travel in nearly circular orbits around 
the galactic center. The determination of the rotation curve 
Vc{iaj), where zu is the galactocentric radius, and in par- 
ticular Vi;{vjq) is one of the hardest problems in Galactic 
Structure (see section 2.3 for the local shape of the rotation 
curve and the determination of Vc{zuq)). The rotation curve 
is determined by observations of the kinematics of the gas, 
and in particular of the neutral hydrogen 21 cm line, but the 
rotation curve is not well established for the galactic radii 
zu > zuq. Observations of outer spiral galaxies indicate that 
the rotation curve remains more or less flat after attaining 
its maximum (e.g. Casertano & van Gorkom 1991): the ro- 
tation curve of a mass model of the Milky Way therefore is 
likely to behave similarly. 



2.3 Oort constants and local circular speed 

As tlie global sliape of the rotation curve is not known 
precisely, the determination of its local shape in the solar 



Vr = Adsm2l 

vt = Adcos2/ + Bd 



with 



A 



2^zu 



dtu^"0' 



B 



1 ,«£ dVc . 

'2^zj ^ dzj'^'^ 



(1) 
(2) 

(3) 



Kuijken & Tremaine (1991) showed that the Taylor ex- 
pansion terms arising from non-circularity of the orbits are 
negligible, just as tliey should be if the Milky Way is ax- 
isymmetric. Oort (1927b) showed for the first time this si- 
nusoidal effect of the galactic rotation on the radial velocities 
and on the proper motions. He found A = 19kms~^kpc~^ 
and B = — 24km s^^kpc"^. Indeed, by fitting Eq.(2) to ob- 
served proper motions, one can determine A and B if one 
is sure that the frame is not rotating. This last require- 
ment was pretty unsure before the Hipparcos mission. The 
most reliable determination of the Oort constants based 
on the proper motions of the Cepheids that were mea- 
sured by Hipparcos has been derived by Feast & Whitelock 
(1997) who found A = 14.82 ± 0.84km s"^kpc"^ and B = 
— 12.37 ± G.64kms~^kpc~^. These values were confirmed by 
Mignard (2000) who found A = 14.5 ± 1.0km s^'kpc"^ and 
B = — 11.5±1.0km s~^kpc~^ by using proper motions of dis- 
tant giants. This indicates that the rotation curve is slightly 
declining in the solar neighbourhood. 

The determination of vdzuQ) follows from the de- 
termination of the Oort constants: using the values of 
Feast & Whitelock (1997), we find Vc{zjq) = (218 ± 
8kms^^)(ti70/8kpc), a value which is consistent with 
t)c(o7o) = 220kms~^ that we choose to adopt in this paper. 
Since we impose Vc{zjq) = 220 kms^^ and zjq = 8±0.5 kpc 
for all the potentials defined in section 3, the value of A—B = 
^'"^"^"^ is in the fixed interval 27.6 ± 1.7km s"^kpc"\ This 
interval is in accordance with the one deduced from the val- 
ues of A and B determined by Feast & Whitelock (1997), 
i.e. ^^-^^ = 27.2 ± 0.9km s"^kpc"\ So, the only relevant 
constraint relative to the Oort constants for our potentials 
will be the value of ^(zuq). 

However, The Oort constants are not considered as a 
very strong constraint in this paper. Indeed, the measure- 
ment of the proper motion of the compact radio source 
Sgr A* (Backer 1996), seems to indicate that A — B = 
30.1±0.8km s^'kpc^' when assuming that Sgr A* is station- 
ary with respect to the galactic center. This is inconsistent 
with the determination based on Hipparcos data, and leaves 
us with an uncertainty which still awards a resolution. 
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2.4 Local dynamical mass 

The mass density in the solar neighbourhood pQ is an essen- 
tial constraint for any mass model. It can be surmised that 
this parameter is not directly observed: it is deduced from 
the positions and velocities of tracer stars in tlic direction 
perpendicular to the galactic plane (the ^-direction), using 
the Boltzmann and Poisson equations in various forms. 

Jeans equations (first order moments of the Boltzmann 
equation) imply that the vertical acceleration is related 
to the vertical velocity dispersion and to the vertical num- 
ber density g{^z) of a population of stars by: 

2d/n(5(^)/g(0)) 



(4) 



if the vertical motion can be separated from the radial and 
azimuthal motion of the stars (Oort-Lindblad approxima- 
tion), if the velocity ellipsoid is aligned with the cylindrical 
coordinate axes (i.e {v^Vz) = 0), and if the population is 
isothermal (i.e. al is constant as a function of z). 

Oort (1932; 1960) applied this now classical formula to 
late-type stars (A to M) with the assumption that they were 
old enough to have become dynamically well mixed in the 
z-direction. Then he derived pQ using Poisson equation for 
nearly circular orbits: 



A^Gp^ = -^-2{A'-B') 



(5) 



Oort found p© = O.O9M0pc"' in 1932 and = 
O.ISMqpc^^ in 1960: this last result indicated that there 
might be a lot of dark matter in the disc, by comparison 
with star counts. 

Afterwards, many other similar studies (Yasuda 1961; 
Eelsalu 1961; Woolley & Stewart 1967; Turon Lacarrieu 
1971; Gould & Vandervoort 1972; Jones 1972; Balakirev 
1976; Hill, Hilditch & Barnes 1979) gave discrepant results, 
suffering from inhomogeneities in the data, systematic errors 
due to the use of photometric distances, and undersampling 
near the galactic plane. Because of this undersampling, some 
studies were even based on young O and B stars, assuming 
that the gas and dust out of which they recently formed was 
already roughly relaxed (Stothers & Tech 1964; von Hoerner 
1966). 

More recently, Bahcall (1984a,b,c) described the disc 
matter as a series of isothermal components and analyzed 
the nonlinear self-consistent equations in which the mat- 
ter produces the potential (Poisson equation) and is also 
afi^ected by the potential (Jeans equations for each isother- 
mal component). Bahcall, Flynn & Gould (1992) applied 
this method to a sample of K giants and found pQ = 
0.26Mopc~^, a result leading once again to the presence 
of dark matter in the disc. 

At about the same time, Kuijken & Gilmore (1989a,b,c, 
1991) used another method which does not use the Jeans 
equations, inspired by the method of von Hoerner (1960). 
Tills method is based on the assumption that the phase 
space distribution function Fz{z,Vz) of any tracer popula- 
tion depends only on 



Ez 



Vz{z) + -vi 



where 
dVzjz) 
dz 



K.{z). 



(6) 



(7) 



This assumption follows from the classical separability of the 
vertical motion of the stars and from the Jeans (1915) the- 
orem in one dimension. The density in configuration space 
g{z) of a tracer population is then related to its density in 
phase space by the integral equation: 



9{z) 



f 



Fz{z,Vz)dvz = 2 



f 



Fz{Ez 



V^iEz - Vz) 



dEz 



f{Vz) = f{Vz{z)) 



(8) 



So, there is a unique relation between g{z) and Fz(Ez): Kui- 
jken & Gilmore (1989c, 1991) inverted the Abel transform 
(8) and used the space density profile g(z) of distant K stars 
to predict their velocity distribution at different heights for 
different Vz{z). Then they compared these predictions to the 
velocity data and used a maximum likelihood technique to 
select the best-fitting potential. The data they used were too 
far from the plane to constrain and they found a surface 
mass density between z = ±1.1 kpc of 71M0pc~^, a result 
rejecting the presence of dark matter in the disc. 

So, all the investigations between 1932 and the Hippar- 
cos mission had failed to converge to a reliable determination 
of pq: they were all very uncertain, essentially because of 
inhomogeneities in the tracer samples, undersampling near 
the equatorial plane and the use of photometric distances. 
Hipparcos data solved all these problems: the completeness 
of the Hipparcos survey for stars brighter than m^. = 8.0 
solved the homogeneity problem, the dense probe near the 
plane eliminated the undersampling, while the accurate par- 
allaxes eliminated the use of photometric distances. 

Crczc et al. (1998) and Holmberg & Flynn (2000) used 
Eq. (8) to determine p© from complete samples of nearby 
A-F stars. Given the functions g{z) and Fi{z = 0,v^^o) = 
Fz{Ez) from the observed vertical density and velocity dis- 
tribution, the function Vz{z) can be derived. Then, to deter- 
mine Pq, the Oort constants have to be used in the Poisson 
equation (5): these are also much better known since the 
Hipparcos mission (see section 2.3). Creze et al. (1998) esti- 
mated po = 0.076 ±O.O15M0pc"^ while Holmberg & Flynn 
(2000) estimated po = 0.102 ± O.OlOMopc"^ These dif- 
ferences between local density estimates using almost the 
same Hipparcos data are due to different a priori hypothesis 
on the vertical potential Vz - Recent investigations at Stras- 
bourg University using a local Hipparcos sample combined 
with two samples at the galactic poles confirm the values 
obtained by Holmberg & Flynn (2000) (Siebert, Bienayme 
& Soubiran 2002). 

We conclude that any mass model of the Milky Way 
must have its solar neighbourhood mass density in the range 
O.O6M0pc"^ < P0 < O.12M0pc"^ in order to match Hip- 
parcos latest findings (rejecting the presence of a disc-like 
dark matter component). However, it should be remarked 
that all the results presented here are based on the Oort- 
Lindblad approximation, which may not be valid: therefore 
we do not exclude that another subset of potentials, that do 
not match our imposed constraints based on Oort-Lindblad 
approximation, could also be useful for dynamical modeling 
of the Milky Way. 
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3 THE POTENTIALS 

The Stackel potentials are non-rotating potentials for which 
the Hamilton-Jacobi equation is separable and for which all 

orbits admit three analytic integrals of the motion (Stackel 
1890): they form the most general set of potentials that con- 
tain one free function, for which three exact integrals of the 
motion are known, and which can be relevant as models for a 
global potential in galactic dynamics (Lynden-Bell 1962). In 
this section, we present the mathematical form of the class 
of Stackel potentials that are studied in this paper. 



3.1 Coordinate system 

Axisymmetric Stackel potentials are best expressed in 
spheroidal coordinates {X,(f>,u), with A and u the roots for 
r of 



+ 



1 



Q < 7 < 0, 



(9) 



T + a T + "/ 

and (c7, (/), z) cylindrical coordinates. The parameters a and 
7 arc both constant and we assume them smaller than zero. 
It is convenient to define the axis ratio of the coordinate 
surfaces as e = ^ with a = —a^ and 7 = — c^. Together with 
the focal distance A = ^7 — a, the axis ratio defines the 
coordinate system. 



3.2 Three-component Stackel potentials 

An axisymmetric potential is of Stackel form if there exists a 
spheroidal coordinate system (A, (p, v) in which the potential 
can be written as 



/(A) - f{y) 



(10) 



for an arbitrary function /(t) = {t + 'y)G{T), G > 0, r = 
A, V. The function —G{X) then represents the potential in 
the z = plane. 

The Milky Way is composed of several mass compo- 
nents: the bulge, the thin disc, the thick disc, the stellar 
halo, the dark halo, and the (dynamically insignificant) in- 
terstellar medium. It is not fundamental that a mass model 
aknowledges explicitly the existence of each of these compo- 
nents. For example, BD presented a set of two-component 
(halo-disc) Stackel potentials with a fiat rotation curve. 

Our goal is to show that different Stackel potentials are 
able to fit the latest estimates for the fundamental parame- 
ters of the Galaxy. We first generalize the potentials of BD 
by adding a thick disc to them since its existence as a sep- 
arate stellar component is now well established (Ojha et al. 
1994; Chen et al. 2001). Our potentials have thus three mass 
components: two "flat" components and one spheroidal. The 
spheroidal component accounts for the stellar and dark halo, 
and we shall see in section 4 that our potentials turn out to 
have an effective bulge, which enables us to avoid the explicit 
introduction of a bulge-component. 

We assume that all three components of our potentials 
generate a Stackel potential, with three different coordinate 
systems but the same focal distance. It is straightforward to 
show that the superposition of three Stackel potentials is still 
a Stackel potential when all three coordinate systems have 
the same focal distance. Although the functions /thin, /thick 



and /halo are arbitrary, we assume that they each generate 
a Kuzmin-Kutuzov potential, defined by 

GM 



G(r) 



Vt + c 



(11) 



with M the total mass of the system. Such a potential be- 
comes a point mass potential {V = —^^) in the Galactic 
Plane when A ^ 00. We use this potential essentially be- 
cause it is an extremely simple but representative Stackel 
potential: we will show that it is not necessary to use com- 
plicated Stackel potentials in order to match all the known 
and most recently determined parameters of the Milky Way. 

Near the center, in a meridional plane, the lines of con- 
stant mass density corresponding to a Stackel potential are 
approximately ellipsoidal (de Zeeuw, Peletier & Pranx 1986). 
Por a Kuzmin-Kutuzov potential (see e.g. Dejonghe & de 
Zeeuw 1988), when a > c, the isodensity surfaces are flat- 
tened oblate spheroids, and increasing e = 7 produces more 
flattening. So, the ratio e has to be high for the thin disc, 
intermediate for the thick disc and close to unity for the 
halo. 

We first define a class of dimensionless potentials Vp in 
dimensionless units {zUp,Zp), with a focal distance A = 1 
for all three coordinate systems and the central value of 
the potentials equal to —1. Each of these potentials is a 
superposition of three Kuzmin-Kutuzov potentials in three 
different coordinate systems: 

Vp(\thin, Athick, Ahalo- 'Aliin. 'Ahick, I'halo) = 

, /thin(Ailiin) — /lliin ('Aliin) 

"-thin 7 

Athin 'Ahin 

/thick (Athick) — /thick (I'thick) 



-fcthic 



-(1 — /Sthin - /Sthick) 



Athick ~ I'thick 

/halo (Ahalo) — /halo(t^alo) 



(12) 



Ahalo ^^alo 

This new class of potentials is thus defined by five param- 
eters (the three axis ratios of the coordinate surfaces and 
the relative contribution of the thin and thick disc masses 
to the total mass, i.e. fethin and fcthick), which is a reasonable 
augmentation with respect to the BD potentials that were 
defined by three parameters. 
We denote 

athin — athick = 7thin — 7thick = Qi > 

aUiin — ahalo = 7thin — 7halo = ^'2 > ^1 > 0. (13) 

So, we can express the class of potentials Vp as a function 
of Athin, I'thin and the two constants qi and 52 (we also use 
Eq. (11)) to give the final form of Vp-. 



V^(Athin,t'thin) = -GM( 



V Athin + y^'thin 



fethic 



V^thin + qi+ Vthin + qi 
1 ^thin ^thick 



VAthin + 92 + V^thin + 92 



(14) 



The dimensionless rotation curve corresponding to such 
potentials is given by: 

^thin 



GMvOpi 



VAthin(v'A 

thin + Cthin)^ 



© 2002 RAS, MNRAS 000, 1-13 



Three- component Stdckel potentials satisfying recent estimates of Milky Way parameters 5 



A:thic 



(15) 



V-'^thin + qiW^Vhin + qi + \/c?hin + I^Y 

1 ^thin ^thick ^ 
V-'^thin + q2{\/\th\n + 92 + \/c?hin + 92)^ 

where Ti7p denotes the dimensionless galactocentric radius. 
We shall impose constraints on the shape and flatness of the 
rotation curve in section 3.3. 

In order to transform these dimensionless potentials 
into dimensional ones for the Milky Way, we denote the di- 
mensionless radius where the rotation curve attains its first 
maximum as Wp^M- since the Milky Way attains its global 
amplitude of 220 km/s for the first time at a radius of about 
1.5 kpc (Fich & Tremaine 1991), we define a distance scale 
factor rs = The conversion between dimensionless 

and dimensional distances is then given by: 



tx7(kpc) 
2;(kpc) 



rsZp 



(16) 



Then, the total mass M of the Galaxy is adjusted in 
such a way that the dimensional circular velocity at the 
solar radius {zuq = 8 ± 0.5 kpc, see section 2.1) is equal 
to 220kms~ : we obtain thus a minimum and a maximum 
value for M, for the two extreme values of the galactocen- 
tric distance of the sun. That adjustment also fixes the local 
mass density in the solar neighbourhood pq (see section 3.3). 



3.3 Selection criteria 

We shall now establish the features that a potential (as de- 
fined in section 3.2) must have to be considered as a plausible 
potential for the Milky Way, in the light of the observational 
constraints reviewed in section 2. 

By definition of the dimensional potentials (see section 
3.2), the local circular speed VcivJc-)) is equal to 220kms~^ 
for all the potentials. The first fundamental selection crite- 
rion is the fiatness of the rotation curve: this feature can 
be examined in the dimensionless frame-work. Even the BD 
potentials, with only two mass components, could produce 
many difi^erent shapes of rotation curves. So, we adopt the 
same simple diagnostic as BD: we denote izip^M the dimen- 
sionless radius where the circular speed attains its first max- 
imum and we look for a range in vj where Vc{vj) remains 
larger than 80% of the maximum velocity and is thus more 
or less constant. We denote vjp^F the dimensionless radius 
where vdvop^p) = 0.8t)c (rop^M). A rotation curve is consid- 
ered sufficiently fiat if: 



^p,F — T^p^M 



> 



(17) 



This is a minimum requirement. 

The second selection criterion is based on the latest de- 
terminations of the Oort constants. For all our potentials. 



fc(ro0) 



27.6 ± 1.7km s"^kpc"^ 



(18) 



which is in accordance with the value determined by Feast 
& Whitelock (1997), ""^J^^ = 27.2 ± O.gkms^ikpc^i. The 
first derivative of the circular velocity in the solar neighbour- 
hood corresponding to the Oort constants found by Feast & 
Whitelock (1997) is ^(073) = -2.4±1.2kms"ikpc"^ Our 



potentials have two extreme values for this derivative, de- 
pending on the position of the sun; in order to fit the above 
interval, we select the potentials such that: 



max(^(€a0)) > -3.6kms-^kpc-^ 

min(^(ti7(5)) < -1.2km s-^kpc"^ 



(19) 



This feature is not as essential as the flatness of the rotation 
curve, because of the intriguing measurement of the proper 
motion of Sgr A* (see section 2.3). 

The last selection criterion is the local mass density in 
the solar neighbourhood: this number is determined by the 
adopted total mass M. We look for its values in both ex- 
treme positions for the sun {tuqjZq) = (7.5, 0.004) kpc and 
{■ojqjZq) = (8.5, 0.02) kpc. Following the Hipparcos latest 
findings reviewed in section 2.4, we select the potentials such 
that: 



max(p0) > O.O6M0PC 
min(p0) < O.I2M0PC" 



(20) 



This feature is more important than the Oort constants: we 
give the priority to that constraint, contrarily to Sevenster 
et al. (2000) who studied the structure of the inner Galaxy, 
and therefore constructed Stackel potentials with reasonable 
values for the Oort constants but values for pQ that are too 
low. 

The total mass of the Galaxy, the mass fractions of the 
discs and the fiattening and scale length of the components 
are not well established observationally and are not consid- 
ered as fundamental constraints for the potential. 



4 RESULTS 

Our goal is now to find and select, among the class of po- 
tentials defined in section 3, some representative potentials 
that difl^er with respect to their form and features, and that 
satisfy the fundamental constraints reviewed in section 2 
and quantified in section 3.3. As it is difficult to visualize 
a five-dimensional parameter space, we shall first visualize 
its structure when tthin and ehaio are fixed (the "winding 
staircase"). Then we shall restrict the parameter space by 
imposing additional constraints on the fiatness of the thick 
disc. We shall finally select five representative potentials, 
listed in Table 5. 



4.1 The "winding staircase" 

As an example of the consequences of the choice of the se- 
lection criteria of section 3.3, we look for all the values of 
ethick, fcthin and fethick that yield potentials satisfying the se- 
lection criteria for ethin = 75 and ehaio = 1-02 (with the 
thick disc always thicker than the thin disc, i.e. Cthick < 75). 
The accordance with the selection criteria results from a 
precise mixing of the two discs and of the halo. Figure 1 
illustrates the volume in parameter space corresponding to 
these satisfactory potentials: the volume looks like a "wind- 
ing staircaise". When fethick = 0, we see that all the values 
of ethick are allowed when 13% < fcthin < 15%, which results 
from the fact that no thick disc is in fact present. When 
fcthin is close to its maximum possible value (15%), fcthick 
has to be zero except when ethick is very close to 1, i.e. when 
the thick disc is a pretty round component (similar to the 
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k thick 




0.04 

0.06 

0.08 

k thin O-"" 



40 

a thick/c thick 



0.16 



Figure 1. This "winding staircase" figure displays the possible values of the coordinate axis ratio ^ = ^ for the thick disc and the 
possible values of the contributions k of the discs to the total mass in order to satisfy the selection criteria of section 3.3 (for fixed values 

Ethin = 75 and fhalo = 1-02, and with the thick disc always thicker than the thin disc, i.e. Cthick < 75). It gives a rough vision of the 
region of parameter space that satisfies the criteria. Only the region 1.3 < ethick < 2 is in fact relevant for a model of the Milky Way 
(see section 4.2). 



halo) and helps to keep the rotation curve flat. When fethin 
is smaller than 15%, the possibilities for fcthick are more nu- 
merous, i.e. the thick disc can take a part of the mass. When 
fcthin attains the critical value of 12%, the volume is inflected 
and the thick disc has to be thin and non-zero in order to 
counter-balance the lack of mass in the thin disc. Decreas- 
ing the mass of the thin disc forces the thick disc to become 
thinner and more massive, yielding the "winding staircase" 
volume of Figure 1. The similar volumes in parameter space 
for ethin > 75 have the same form, and are bigger essentially 
because there is more freedom for ethick. 



4.2 Constraints on the scale height of the thick 
disc 

In Figure 1, there arc some solutions with a thick disc more 
massive than the thin disc. Even though the mass fractions 
of the discs and the flattening of the components are not 
well established and should be tested in a dynamical study, 
we know that the mass fraction of the thick disc is smaller 
than that of the thin disc (and represent at most 13% of the 
local thin disc density in the solar neighbourhood) and the 
latest determination of the thick disc scale height based on 
star count data from the Sloan Digital Sky Survey is 665pc 
(Chen et al. 2001). However, Chen et al. (2001) insist on 



the difficulty to converge to a definitive answer: some other 
studies indicate that the scale height could be of the order of 
1 kpc (Gilmore 1984; Ojha et al. 1996). We shall reject the 
potentials that are completely inconsistent with those char- 
acteristics (0.6 kpc < hz < 1 kpc) and in particular those 
with a thick disc more massive than the thin disc. 

In order to determine the interval of axis ratios that 
should be considered for the thick disc, we have fitted 
the vertical mass distribution corresponding to a simple 
Kuzmin-Kutuzov potential with rg = 1 to an exponential 
law e~^^^'. We conclude that the potentials with 1.3 < 
ethick < 2 have a scale height between 665 pc and 1 kpc 
and are the ones we should examine in detail. For each axis 
ratio. Table 1 gives the corresponding scale height. How- 
ever, Figure 2 reveals that the exponential fit is not valid 
for the axis ratios used to model the thin disc, which is not 
surprising since a thin disc could be better understood as 
a superposition of isothermal sheets than by a simple expo- 
nential law. 



4.3 The final selection 

In this section, we look for some three-component potentials 
with different forms and features, all satisfying the selection 
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2 3 
z(kpc) 





z(kpc) 



z(kpo) 



Figure 2. Profile of the logarithm of vertical density at tjj = 8 kpc for Kuzmin-Kutuzov potentials with = 1 and e — 1.3 (top left), 
e = 2 (top right), e = 75 (bottom left) and e = 200 (bottom right). Only the two first cases resemble exponentials. 



criteria defined in section 3.3 and consistent with what is 
known about the thick disc. 

First of all, we look for the two-component BD po- 
tentials that satisfy the new selection criteria: they are 
listed in Table 2. All the two-component potentials with 
CdiBc = 50, 75, 130, 200 and ehaio = 1.005, 1.01, 1.02, 1.03 are 
examined. We do not consider discs with a/c > 200 because 
then the uncertainty on pQ becomes too large. We see in 
Table 2 that, in order to reproduce the Oort constants in 



the two-component framework, the shape of the halo can- 
not vary (ehaio = 1.02). 

If we take the two-component potentials as a starting 
point, there are two ways to add a thick disc. The first way 
is to decrease the contribution of the halo and to put the 
remaining mass into the thick disc: the local density in the 
solar neighbourhood is then slightly larger while the rota- 
tion curve is decreasing faster. If we take the third potential 
of Table 2 as a starting point, the first potential of Table 
5 (potential I) illustrates this first case. The other way is 



© 2002 RAS, MNRAS 000, 1-13 



8 B.Famaey & H.Dejonghe 




Figure 3. Mass isodensity curves in a meridional plane for the five potentials of Table 5, with two different scales (left panel: zoom on 
the disc, right panel: large scale view of the halo). 
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Table 2. The characteristics of the two-component BD potentials witli ejig^ = 50,75,130,200 and Chalo ~ 1-005,1.01,1.02,1.03 have 
been examined. The potentials satisfying the selection criteria are listed in this table (with a step of 0.01 for the relative contribution of 
the disc). Columns 1 and 2 contain the axis ratios of the coordinate surfaces for the disc and the halo. Column 3 contains the relative 
contribution of the disc to the total mass. Column 4 contains the extent of the flat part of the rotation curve (see Eq. 17). Column 5 
contains the minimum and maximum local spatial density, while column 6 contains the minimal and maximal local radial derivative of 
the circular velocity, each time for the two extreme positions of the sun. 



f-disc Chaio fcdisc Ef P0 in Mqpc 3 ^(ro0)inkms ^kpc 



75 


1 


02 





11 


13.07 





04, 





06 


-1.90, 


-1.18 


75 


1 


02 





12 


11.74 





04, 





06 


-2.03, 


-1.42 


75 


1 


02 





13 


10.43 





04, 





07 


-2.25, 


-1.74 


75 


1 


02 





14 


9.25 





04, 





07 


-2.46, 


-2.05 


75 


1 


02 





15 


8.20 





05, 





08 


-2.68, 


-2.37 


130 


1 


02 





08 


17.16 





04, 





06 


-2.07, 


-1.12 


130 


1 


02 





09 


16.03 





04, 





07 


-1.75, 


-0.85 


130 


1 


02 





10 


14.56 





05, 





08 


-1.73, 


-0.92 


130 


1 


02 





11 


13.03 





05, 





09 


-1.86, 


-1.15 


130 


1 


02 





12 


11.70 





06, 





10 


-2.00, 


-1.39 


130 


1 


02 





13 


10.38 





06, 





11 


-2.22, 


-1.72 


130 


1 


02 





14 


9.20 





07, 





11 


-2.44, 


-2.04 


130 


1 


02 





15 


8.08 





07, 





12 


-2.69, 


-2.38 


200 


1 


02 





08 


17.21 





04, 





09 


-2.02, 


-1.07 


200 


1 


02 





09 


16.07 





05, 





10 


-1.71, 


-0.81 


200 


1 


02 





10 


14.59 





06, 





12 


-1.69, 


-0.88 


200 


1 


02 





11 


13.05 





07, 





13 


-1.82, 


-1.12 


200 


1 


02 





12 


11.64 





07, 





14 


-2.00, 


-1.40 


200 


1 


02 





13 


10.32 





08, 





16 


-2.22, 


-1.72 


200 


1 


02 





14 


9.14 





08, 





17 


-2.44, 


-2.04 


200 


1 


02 





15 


8.08 





09, 





18 


-2.66, 


-2.36 



Table 1. Column 1 contains the axis ratio of the coordinate sur- 
faces for the thick disc. Column 2 gives the corresponding scale 
height for a Kuzmin-Kutuzov potential with rg — 1. 



e 


/iz(pc) 


1.3 


914 


1.4 


834 


1.5 


782 


1.8 


696 


2 


663 



to decrease the contribution of the thin disc and to put the 
remaining mass into the thick disc: the local density is then 
slightly decreasing while the rotation curve is more flat. If 
we take the fourth potential of Table 2 as a starting point, 
the first potential of Table 5 illustrates this second case: this 
potential (potential I) is a three-component potential satis- 
fying the selection criteria and is selected to be analyzed in 
detail and confronted with kinematical surveys in the future. 

The presence of a third component allows more freedom 
for the shape of the halo, so we look for three-component 
potentials with a halo rounder than thaio = 1.02. In order to 
keep the rotation curve flat and retain the local density as 
well as the Oort constants in the allowed interval, we need 
to couple a very thin disc with the rounder halo: indeed, our 
investigations show that no solution can be found for ethin = 
50 and ehaio = 1.01. However, if we take ethin = 200, Table 
3 gives solutions for a halo with ehaio = 1.01: we select the 
solution where the mass of the thick disc relative to the thin 
disc is the smallest (potential III of Table 5). Remark that a 
similar Table for Chaio = 1-02 would contain 292 entries and 



is omitted here. There are much less solutions when ethin = 
75, as can be seen in Table 4. However, as stated in section 
3.3, wc do not assign high priority to the Oort constants, 
and we select a potential with ethin = 75, ehaio = 1-01 and a 
relative mass of the thick disc relative to the thin disc of 13% 
(i.e. a smaller fraction than any of the solutions of Table 4), 
but with a quite large local radial derivative of the circular 
speed (potential IV of Table 5). 

For ehaio = 1.005, it is totally impossible to find a po- 
tential satisfying the Oort constants criterion: the radial 
derivative of the circular speed in the solar neighbourhood 
is always positive. Nevertheless, if one is willing to ignore 
the estimates of A and B, one could select a potential with 
ehaio = 1.005 and reasonably low values for the cirular speed 
radial derivative in the solar neighbourhood (potential V of 
Table 5). 

Finally, we select a potential (potential II) satisfying 
all the criteria, for which the interval in p© is precisely 
[O.O6Mopc"'',O.12M0pc"^], and which is close to the Chen 
et al. (2001) findings , i.e. ethin = 200 which is the thinnest 
thin disc that we consider in order not to have a too large 
interval for p©, a relative mass of the thick disc to the thin 
disc of 10%, a scale height of the thick disc of 612.5 pc and 
a relatively large Ef (Ef = 13.44). 

Table 5 summarizes the main features of the selected 
Stackel potentials with different forms and features and that 
we shall use for dynamical modeling of the Milky Way: po- 
tentials HI and V have a very thin disc associated with 
a quite massive thick disc and a pretty round halo, while 
potentials I and IV have a thicker thin disc with a quasi- 
negligible thick disc (Figure 3 shows the mass isodensity 
curves of each potential in a meridional plane for two differ- 
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Table 3. The three-component potentials with fthin = 200 and fhalo = 1-01 that satisfy the selection criteria are hsted in this table 
(with a step of 0.01 for the relative contribution of the two discs). Columns 1, 2, 3 contain the axis ratios of the coordinate surfaces for 
the two discs and the halo. Columns 4 and 5 contain the relative contribution of the repectively thin and thick disc to the total mass. 
Column 6 contains the extent of the flat part of the rotation curve (see Eq. 17). Column 7 contains the minimum and maximum local 
spatial density, while column 8 contains the minimal and maximal local radial derivative of the circular velocity, each time for the two 
extreme positions of the sun. 



ethin <;thick Chaio hhin ^'thick Ep PQ in Mqpc 3 ^(ro0)inkms ^kpc 



200 


1 


3 


1 


01 





09 





07 


10.69 





06, 





12 


1 


34, 


-1.13 


200 


1 


3 


1 


01 





09 





08 


10.06 





06, 





12 


-1 


64, 


-1.45 


200 


1 


3 


1 


01 





10 





06 


9.51 





07, 





15 


1 


35, 


-1.29 


200 


1 


3 


1 


01 





10 





07 


8.95 





07, 





14 


-1 


64, 


-1.59 


200 


1 


3 


1 


01 





10 





08 


8.45 





07, 





13 


-1 


92, 


-1.88 


200 


1 


3 


1 


01 





11 





04 


8.51 





09, 





18 


-1 


22, 


-1.17 


200 


1 


3 


1 


01 





11 





05 


8.02 





08, 





17 


-1 


51, 


-1.45 


200 


1 


4 


1 


01 





09 





07 


9.81 





06, 





13 


1 


46, 


-1.32 


200 


1 


4 


1 


01 





09 





08 


9.11 





06, 





12 


-1 


77, 


-1.65 


200 


1 


4 


1 


01 





10 





06 


8.63 





07, 





14 


1 


50, 


-1.47 


200 


1 


4 


1 


01 





10 





07 


8.07 





07, 





14 


-1 


78, 


-1.78 


200 


1 


5 


1 


01 





08 





07 


10.86 





06, 





11 


-1 


25, 


-1.00 


200 


1 


5 


1 


01 





09 





06 


9.93 





07, 





13 


-1 


21, 


-1.09 


200 


1 


5 


1 


01 





09 





07 


9.10 





06, 





13 


-1 


54, 


-1.45 


200 


1 


5 


1 


01 





09 





08 


8.34 





06, 





12 


-1 


86, 


-1.80 


200 


1 


5 


1 


01 





10 





06 


8.01 





07, 





15 


1 


59, 


-1.57 


200 


1 


8 


1 


01 





08 





07 


9.42 





06, 





12 


-1 


34, 


-1.22 


200 


1 


8 


1 


01 





09 





06 


8.55 





07, 





13 


1 


34, 


-1.32 


200 


2 




1 


01 





08 





07 


8.75 





06, 





12 


-1 


38, 


-1.31 



Table 4. The three-component potentials with tthin = 75 and Chaio = 1-01 that satisfy the selection criteria are listed in this table (with 
a step of 0.01 for the relative contribution of the two discs). Columns have the same meaning as in table 3. 



Cthin 


Cthick 


Chalo 




^^thick 




Pq in Mqpc 3 


^(ro©) in kms ^kpc ^ 


75 


1.3 


1.01 


0.11 


0.05 


8.21 


0.05, 0.07 


-1.49, -1.45 


75 


1.4 


1.01 


0.11 


0.04 


8.06 


0.05, 0.07 


-1.32, -1.27 


75 


1.5 


1.01 


0.10 


0.05 


8.87 


0.04, 0.07 


-1.26, -1.23 



ent scales). It should be noted that the total masses asso- 
ciated with those potentials are very different and become 
larger with a rounder halo and that a rounder halo implies 
that this halo is much more extended. A closer look to the 
mass density in the equatorial plane indicates that, for each 
potential, the density grows faster than an exponential in 
the central 3 kpc corresponding to the bulge region: the po- 
tentials have thus an effective bulge, which did enable us to 
avoid the introduction of an explicit bulge component. We 
have fitted the mass density in the plane to an exponen- 
tial law down to = 3 kpc in order to check that the scale 
length of the disc is realistic: the last column of Table 5 gives 
the scale length corresponding to each selected potential and 
we conclude that they are realistic but do not distinguish the 
different potentials. For the potentials with the biggest and 
smallest scale length, Figure 4 illustrates the shape of the 
logarithm of the density in the plane (the other potentials 
have a similar shape for that curve). Finally, Figure 5 shows 
the rotation curve associated with each of the five selected 
potentials: the rotation curve of potential II is more flat than 
the one of potential I in the vicinity of the sun, while the ro- 
tation curves of the potentials with Chaio = 1.01 (potentials 
III and IV) are even more fiat and the one of potential V is 
slightly increasing. 



5 CONCLUSIONS 

In this paper, we have shown that some different simple 
Stackel potentials can fit most known parameters of the 
Milky Way (especially Hipparcos latest findings). First, we 
have reviewed the galactic fundamental parameters that 
any Milky Way potential must match and that investiga- 
tions following the Hipparcos mission have readjusted. Then, 
we have generalized the two-component BD potentials by 
adding a thick disc, and we have studied how the parame- 
ters can vary in order to satisfy selection criteria based on 
the latest observational constraints. We have shown that the 
presence of a thick disc allows more fiexibility in the form 
of the potentials, especially for the shape of the halo and 
we have selected five different valid potentials listed in Ta- 
ble 5. It should be noted that, in fact, there could be two 
different thick discs in the Galaxy, a very thick one (Chiba 

6 Beers 2000; Gilmore, Wyse & Norris 2002) and a smaller 
one (Soubiran, Bienayme & Siebert 2002): in that case, the 
three-component modelling presented in this paper could be 
easily extended, but this would imply a growth of parameter 
space. 

A major result of this paper is that, even though Stackel 
potentials are negligible in the set of all potentials, many of 
them are still able to match the most recent estimates for the 
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Table 5. Among the class of potentials defined in section 3, five different potentials regarding form and features have been selected. 
Columns 1, 2, 3 contain the axis ratios of the coordinate surfaces for the two discs and the halo. Columns 4 and 5 contain the relative 
contribution of the repectively thin and thick disc to the total mass. Column 6 contains the extent of the fiat part of the rotation curve 
(see Eq. 17). Column 7 contains the scale factor which corresponds to the focal distance of the coordinate system of the dimensional 
potential. Column 8 contains the minimum and maximum local spatial density in Mqpc^^, while column 9 contains the minimum and 
maximum local radial derivative of the circular velocity in kms~^kpc~^ and column 10 the minimum and maximum total mass of the 
Galaxy in IO^^Mq, each time for the two extreme positions of the sun. The scale length in the equatorial plane down to 3 kpc has been 
calculated and is presented in column 11 (in kpc). 









fhalo 






Ef 








M 


scale length 


I 


75 


1.5 


1.02 


0.13 


0.01 


9.86 


0.93 


0.04, 0.07 


-2.51, -2.05 


2.37, 2.41 


2.73 


II 


200 


1.8 


1.02 


0.10 


0.01 


13.44 


0.88 


0.06, 0.12 


-2.05, -1.31 


2.37, 2.41 


2.63 


III 


200 


1.3 


1.01 


0.11 


0.04 


8.51 


0.95 


0.09, 0.18 


-1.22, -1.17 


3.19, 3.22 


2.65 


IV 


75 


1.8 


1.01 


0.11 


0.015 


9.07 


0.98 


0.05, 0.08 


-0.62, -0.55 


3.56, 3.58 


2.78 


V 


200 


1.3 


1.005 


0.07 


0.01 


18.30 


1.01 


0.11, 0.23 


+0.69, +1.47 


6.13, 6.20 


2.72 



T 




galactic radius(kpc) 



Figure 4. The logarithm of the mass density in the equatorial plane for the two potentials with extreme scale lengths. These curves very 
much resemble each other, and the effective bulge appears clearly. 
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Figure 5. The rotation curves of the five selected potentials of Table 5. The total mass used to plot thes curves is the mean total mass 
of the two extreme values of Table 5. 



parameters of the Milky Way, and furthermore very simple 
ones (superpositions of three Kuzmin-Kutuzov potentials) 
are sufficient to do this. 

The potentials defined in this paper have already been 
used in Famaey et al. (2002). Wc plan to further validate 
each of the five proposed potentials by confronting them 
with kinematical surveys. This we shall do by constructing 
three-integral analytic distribution functions in those poten- 
tials, using a quadratic programming technique (Dejonghe 
1989; for an overview see Dejonghe et al. 2001). 
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